function [CSPepoch,bci]=bciCSPtransform(bci,epoch,label)

if isempty(bci.param.CSPbandFreqs),
    bandFreqs=[0 0];    
else
    bandFreqs = bci.param.CSPbandFreqs;
end
    
if nargin==3, % csp estimation,
    % determine classes and combinations
    classes = unique(label);
    nClasses = length(classes);
    nPairs = (nClasses-1)*nClasses/2;
    classPair=zeros(nPairs,2);
    
    k=0;
    for k1=1:nClasses-1,
        for k2=k1+1:nClasses,
            k=k+1;
            classPair(k,:) = [k1,k2];
        end
    end
    
    bci.cspW=zeros(size(epoch,1),size(epoch,1),size(bandFreqs,1),nPairs);
    CSPepoch = zeros(size(epoch,3),bci.param.CSPfilter*2,size(bandFreqs,1),nPairs);
    cspOvsOmask=zeros(bci.param.CSPfilter*2,size(bandFreqs,1),nPairs);
    filtEpoch = zeros(size(epoch));

    % estimate csp for each band and each class pair
    for band = 1:size(bandFreqs,1),
        % filter data
        [B,A]=bciCSPtransform_getBA(bandFreqs(band,:),bci.init.filtSrate);
        if ~isempty(B),
            for tr=1:size(epoch,3),
                filtEpoch(:,:,tr) = filtfilt(B,A,epoch(:,:,tr)')';
            end
        else
            filtEpoch = epoch;
        end
        for c = 1:nPairs,
            idx=label==classes(classPair(c,1))|label==classes(classPair(k,2));
            %epochMean = repmat(mean(filtEpoch(:,:,idx),2),[1,size(filtEpoch,2),1]);
            bci.cspW(:,:,band,c) = bciCSPEstimate(filtEpoch(:,:,idx),...-epochMean,...
                                                  label(idx));
            CSPepoch(idx,:,band,c)=transpose(...
                                bciCSPApply(filtEpoch(:,:,idx),...-epochMean,...
                                bci.cspW(:,:,band,c),bci.param.CSPfilter));
            cspOvsOmask (:,band,c) = c*ones(1,size(CSPepoch,2));
        end
    end
    CSPepoch = reshape(CSPepoch,[size(CSPepoch,1),size(CSPepoch,2)*size(CSPepoch,3)*size(CSPepoch,4)])';
    bci.cspOvsOmask = reshape(cspOvsOmask,[1,numel(cspOvsOmask)]);
else % csp application
    nPairs=size(bci.cspW,4);
    CSPepoch = zeros(bci.param.CSPfilter*2,size(bandFreqs,1),nPairs);
    for band = 1:size(bandFreqs,1),
        % filter data
        [B,A]=bciCSPtransform_getBA(bandFreqs(band,:),bci.init.filtSrate);
        if ~isempty(B),
            filtEpoch = filtfilt(B,A,epoch')';
        else
            filtEpoch = epoch;
        end
        for c = 1:nPairs,
            %epochMean = repmat(mean(filtEpoch,2),[1,size(filtEpoch,2),1]);
            CSPepoch(:,band,c)=transpose(...
                                bciCSPApply(filtEpoch,...-epochMean,...
                                bci.cspW(:,:,band,c),bci.param.CSPfilter));
        end
    end
    CSPepoch = reshape(CSPepoch,[1,size(CSPepoch,1)*size(CSPepoch,2)*size(CSPepoch,3)]);
end


function [B,A]=bciCSPtransform_getBA(bandFreqs,srate)

Wn_high = bandFreqs(1)/(srate/2);             % highpass Cutoff Frequenz
Wn_low = bandFreqs(2)/(srate/2);
if Wn_high~=0 && Wn_low~=0, % bandpass
    [B,A] = butter(4,[Wn_high Wn_low]);
elseif Wn_high~=0 && Wn_low==0, % highpass
    [B,A] = butter(4,Wn_high,'high');
elseif Wn_high==0 && Wn_low~=0, % lowpass
    [B,A] = butter(4, Wn_low,'low');
else
    B=[];
    A=[];
end
    
    
    
    
